//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include "ephi.hpp"

#define PRECISION 10000

struct rb : RingBase
{
  rb (const vect3d& pos, prec_t radius, const vect3d& normal) : RingBase(0, pos, radius, 0, normal) { }

  void pos(prec_t t, vect3d& out) { get_pos(t, out); }
  void slope(prec_t t, vect3d& out) { get_slope(t, out); }
};

static void usage (const char *name)
{
  fprintf (stderr, "%s: %s [-r radius] [-s spacing] [-i current] [-o offset]\n", name, name);
  exit(-1);
}

int main (int argc, char *argv[])
{
  prec_t radius = 0.15;
  prec_t spacing = 0.08;
  prec_t current = 1.0;
  prec_t offset = 0.0;

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' || strlen (argv[i]) != 2 )
          usage(argv[0]);

      switch (argv[i][1])
        {
          case 'r': {
            if ( ++i == argc )
                usage(argv[0]);
            radius = atof (argv[i]);
            break;
          }
          case 's': {
            if ( ++i == argc )
                usage(argv[0]);
            spacing = atof (argv[i]);
            break;
          }
          case 'i': {
            if ( ++i == argc )
                usage(argv[0]);
            current = atof (argv[i]);
            break;
          }
          case 'o': {
            if ( ++i == argc )
                usage(argv[0]);
            offset = atof (argv[i]);
            break;
          }
          default: {
            usage(argv[0]);
          }
        }
    }

  printf ("radius = %f m\n", prec2double (radius));
  printf ("spacing = %f m\n", prec2double (spacing));
  printf ("current = %f A\n", prec2double (current));
  printf ("offset = %f degrees\n", prec2double (offset));

  Statics statics;

  prec_t real_r = radius + spacing / PREC_SQRT2;

  vect3d mid(0, 0, 0);
  vect3d cpos;

  ShapedCoilFactory cf;
  cpos = vect3d(0, 0, -real_r); cf.make (statics, cpos, mid - cpos, radius, 0, 1.0, 0);
  cpos = vect3d(-real_r, 0, 0); cf.make (statics, cpos, mid - cpos, radius, 0, 1.0, 0);
  cpos = vect3d(0, real_r, 0); cf.make (statics, cpos, mid - cpos, radius, 0, 1.0, 0);
  cpos = vect3d(real_r, 0, 0); cf.make (statics, cpos, mid - cpos, radius, 0, 1.0, 0);
  cpos = vect3d(0, -real_r, 0); cf.make (statics, cpos, mid - cpos, radius, 0, 1.0, 0);

  // we measure at the 6th coil
  cpos = vect3d(0, 0, real_r);
  vect3d normal(sin((180.0 - offset) * 2 * PREC_PI / 360), 0, cos((180.0 - offset) * 2 * PREC_PI / 360));
  if ( offset == 0.0 ) normal = vect3d(0, 0, -1);
  normal.print("normal");
  rb ring(cpos, radius, normal);

  // length of one segment
  prec_t len = 2.0 * PREC_PI * radius / (prec_t)PRECISION;

  vect3d force(0, 0, 0);
  vect3d torque(0, 0, 0);

  for ( size_t i = 0; i < PRECISION; i++ )
    {
      vect3d pos, slope;
      prec_t t = (prec_t)i/(prec_t)PRECISION;
      ring.pos (t, pos);
      ring.slope (t, slope);
      slope.normalize();
      slope *= len;
      vect3d bf, _ef;
      statics.getFields (pos, bf, _ef);
      vect3d df = slope % bf;
      force += df;
      torque += (pos - cpos) % df;
    }

  force *= current * current;
  torque *= current * current;

  force.print("F");
  torque.print("T");
}

